/* External definitions for single-server queueing system. */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>  /* Header file for random-number generator. */

#define Q_LIMIT 100  /* Limit on queue length. */
#define BUSY      1  /* Mnemonics for server's being busy */
#define IDLE      0  /* and idle. */
#define a (double) 16807.0
#define m (double) 2147483647.0
#define q (double) 127773.0
#define r (float) 2836.0

int next_event_type;
float area_num_in_q, threshhold, area_server_status, mean_interarrival, mean_service, area_lambda, num_delays_required, num_events,
num_in_q, server_status, job_gt_45, sim_time, time_arrival[Q_LIMIT + 1], time_last_arrival, t_star, t_star_last_arrival, time_last_event, time_next_event[3],
total_of_delays, my_random, time_required, time_q_gt1, max_arrival, arrival1, arrival2, currArrival, hour_coming, num_custs_delayed;
FILE  *infile, *outfile;
double hi, test, lo, *seed1, *seed2z;
char arrival, type,terminating_condition;

void  initialize(void);
void  timing(void);
void  arrive(void);
void  depart(void);
void  report(void);
void  update_time_avg_stats(void);
float expon(float mean, double*);
float thinning(void);
float inverse_transform(void);
float randomvar(double*);


main()  /* Main function. */
{
	/* Specify the number of events for the timing function. */

	outfile = fopen("md1.out", "w");
	num_events = 2;

	/* Specify the starting Seed */
	do
	{
		printf("Enter if arrival times arrive with either the (T)hinning or (I)verse transform methods\n");
		scanf("%c", &arrival);
		if(arrival != 'T' && arrival != 'I')
		{
			printf("Incorrect input, try again\n");
		}
	}
	while(arrival != 'T' && arrival != 'I');

	int done = 0;
	seed1 = malloc(sizeof(double));
	seed2 = malloc(sizeof(double));


	do
	{
		printf("Enter terminating condition time(t) or people(p)\n");
		scanf("\n%c", &terminating_condition);
		if(terminating_condition == 't')
		{
			printf("Enter terminating time in min\n");
			scanf("%f", &time_required);
			done = 1;
		}
		else if(terminating_condition == 'p')
		{
			printf("Enter terminating number of people\n");
			scanf("%f", &num_delays_required);
			time_required = -1;
			done = 1;
		}
		else
		{
			printf("Incorrect format try again");
		}
	}
	while(done != 1);
	/* Read input parameters. */

	mean_interarrival = 5;
	mean_service = 4;
	/* Write report heading and input parameters. */
	if(arrival == 'T')
	{
		fprintf(outfile, "M/D/1 System using Thinning method for arrivals with finish requirement specified below\n\n");
	}
	else
	{
		fprintf(outfile, "M/D/1 System using Inverse Transform method for arrivals with finish requirement specified below\n\n");
	}
	fprintf(outfile, "Single-server queueing system\n\n");
	fprintf(outfile, "Mean interarrival time%11.3f minutes\n\n",
			mean_interarrival);
	fprintf(outfile, "Mean service time%16.3f minutes\n\n", mean_service);
	if(terminating_condition == 'p')
	{
		fprintf(outfile, "Number of customers required to finish%14.1f\n\n", num_delays_required);
	}
	else
	{
		fprintf(outfile, "Time required to finish %17.3f\n\n", time_required);
	}
	fflush(stdout);
	/* Initialize the simulation. */
	initialize();

	/* Run the simulation while more delays are still needed. */
	if(terminating_condition == 'p')
	{
		while (1) {
			/* Determine the next event. */
			timing();
			/* Update time-average statistical accumulators. */

			update_time_avg_stats();
			/* Invoke the appropriate event function. */
			if(num_custs_delayed < num_delays_required || next_event_type == 1)
			{
				switch (next_event_type) {
				case 1:
					arrive();
					break;
				case 2:
					depart();
					break;
				}
			}
			else
			{
				break;
			}
		}
	}
	else
	{
		while (1) {
			/* Determine the next event. */

			timing();
			if(sim_time > time_required)
			{
				sim_time = time_required;
			}
			update_time_avg_stats();
			/* Update time-average statistical accumulators. */
			if(sim_time < time_required) {
				/* Invoke the appropriate event function. */
				switch (next_event_type) {
				case 1:
					arrive();
					break;
				case 2:
					depart();
					break;
				}
			}
			else
			{
				num_custs_delayed--;
				break;
			}
			printf("Sim time = %f, Event = %d, total_delay = %f, num_in_q = %f, num_delayed = %f\n", sim_time, next_event_type, total_of_delays, num_in_q, num_custs_delayed);

		}
	}

	/* Invoke the report generator and end the simulation. */

	report();

	fclose(outfile);

	return 0;
}


void initialize(void)  /* Initialization function. */
{
	/* Initialize the simulation clock. */
	int i;
	sim_time = 0.0;

	/* Initialize the state variables. */
	server_status   = IDLE;
	num_in_q        = 0.0;
	time_last_event = 0.0;
	hour_coming = 1.0;
	t_star = 0.0;
	t_star_last_arrival = 0.0;
	/* Initialize the statistical counters. */
	job_gt_45          = 0.0;
	num_custs_delayed  = 0.0;
	total_of_delays    = 0.0;
	area_num_in_q      = 0.0;
	area_server_status = 0.0;
	time_q_gt1		   = 0.0;
	time_last_arrival  = 0.0;
	area_lambda        = 0.0;
	/* Initialize event list.  Since no customers are present, the departure
       (service completion) event is eliminated from consideration. */
	*seed1 = 2062615503;
	*seed2 = 1383670351;
	arrival1 = 10.0/60.0;
	arrival2 = 15.0/60.0;
	currArrival = arrival1;
	threshhold = arrival1*60;
	if(arrival1 > arrival2)
	{
		max_arrival = arrival1;
	}
	else
	{
		max_arrival = arrival2;
	}
	if(arrival == 'T')
	{
		time_next_event[1] = thinning();
	}
	else
	{
		time_next_event[1] = sim_time + inverse_transform();
	}
	time_next_event[2] = 1.0e+30;
}


void timing(void)  /* Timing function. */
{
	int   i;
	float min_time_next_event = 1.0e+29;

	next_event_type = 0;

	/* Determine the event type of the next event to occur. */

	for (i = 1; i <= num_events; ++i)
		if (time_next_event[i] < min_time_next_event) {
			min_time_next_event = time_next_event[i];
			next_event_type     = i;
		}
	/* Check to see whether the event list is empty. */

	/*else
	{
		area_lambda += (min_time_next_event-sim_time) * currArrival;
	}*/

	if (next_event_type == 0) {

		/* The event list is empty, so stop the simulation. */

		fprintf(outfile, "\nEvent list empty at time %f", sim_time);
		exit(1);
	}

	/* The event list is not empty, so advance the simulation clock. */
	//area_lambda += currArrival*(min_time_next_event-sim_time);
	sim_time = min_time_next_event;
}


void arrive(void)  /* Arrival event function. */
{
	float delay;

	/* Schedule next arrival. */
	if(arrival == 'T')
	{
		time_next_event[1] =  thinning();
	}
	else
	{
		time_next_event[1] = inverse_transform();
		time_last_arrival = time_next_event[1];
	}

	/* Check to see whether server is busy. */

	if (server_status == BUSY) {

		/* Server is busy, so increment number of customers in queue. */

		++num_in_q;
		/* Check to see whether an overflow condition exists. */

		if (num_in_q > Q_LIMIT) {

			/* The queue has overflowed, so stop the simulation. */

			fprintf(outfile, "\nOverflow of the array time_arrival at");
			fprintf(outfile, " time %f", sim_time);
			exit(2);
		}



		/* There is still room in the queue, so store the time of arrival of the
           arriving customer at the (new) end of time_arrival. */

		time_arrival[(int)num_in_q] = sim_time;
	}

	else {

		/* Server is idle, so arriving customer has a delay of zero.  (The
           following two statements are for program clarity and do not affect
           the results of the simulation.) */

		delay            = 0.0;
		total_of_delays += delay;

		/* Increment the number of customers delayed, and make server busy. */

		++num_custs_delayed;
		server_status = BUSY;

		/* Schedule a departure (service completion). */

		time_next_event[2] = sim_time + mean_service;
	}
}


void depart(void)  /* Departure event function. */
{
	int   i;
	float delay;

	/* Check to see whether the queue is empty. */

	if (num_in_q == 0) {

		/* The queue is empty so make the server idle and eliminate the
           departure (service completion) event from consideration. */

		server_status      = IDLE;
		time_next_event[2] = 1.0e+30;
	}

	else {

		/* The queue is nonempty, so decrement the number of customers in
           queue. */

		--num_in_q;

		/* Compute the delay of the customer who is beginning service and update
           the total delay accumulator. */
		if(terminating_condition == 't')
		{
			if(sim_time + mean_service < time_required)
			{
				delay            = sim_time - time_arrival[1];
				total_of_delays += delay;
			}
		}
		else
		{
			delay            = sim_time - time_arrival[1];
			total_of_delays += delay;
		}

		/* Increment the number of customers delayed, and schedule departure. */

		++num_custs_delayed;
		time_next_event[2] = sim_time + mean_service;
		if(time_next_event[2] - sim_time + delay > 4.5)
		{
			if(time_required > time_next_event[2])
			{
				job_gt_45++;
			}
			else if(time_required < 0)
			{
				job_gt_45++;
			}
		}
		/* Move each customer in queue (if any) up one place. */

		for (i = 1; i <= num_in_q; ++i)
			time_arrival[i] = time_arrival[i + 1];
	}
}


void report(void)  /* Report generator function. */
{
	/* Compute and write estimates of desired measures of performance. */
	float server_util = area_server_status/sim_time;
	float avg_in_system = (area_server_status+area_num_in_q)/sim_time;
	fprintf(outfile, "Number of customers served%7.1f\n\n", num_custs_delayed);
	fprintf(outfile, "Average delay in queue%11.3f minutes\n\n", total_of_delays / num_custs_delayed);
	fprintf(outfile, "Average number in system%10.3f\n\n", avg_in_system);
	fprintf(outfile, "Average number in queue%10.3f\n\n", area_num_in_q / sim_time);
	fprintf(outfile, "Fraction of Jobs that spent more than 4.5 minutes in the system%10.3f\n\n", (float)job_gt_45/(float)num_custs_delayed);
	fprintf(outfile, "Server utilization%15.3f\n\n", area_server_status / sim_time);
	fprintf(outfile, "Fraction of time that there is more than one customer in queue%10.3f\n\n", time_q_gt1/sim_time);
	fprintf(outfile, "Time simulation ended%12.3f minutes\n\n", sim_time);

}


void update_time_avg_stats(void)  /* Update area accumulators for time-average
                                     statistics. */
{
	float time_since_last_event;

	/* Compute time since last event, and update last-event-time marker. */

	time_since_last_event = sim_time - time_last_event;
	time_last_event       = sim_time;

	/* Update area under number-in-queue function. */

	area_num_in_q      += num_in_q * time_since_last_event;
	if(num_in_q > 1.0)
	{
		time_q_gt1 += time_since_last_event;
	}
	/* Update area under server-busy indicator function. */

	area_server_status += server_status * time_since_last_event;

}


float expon(float mean, double *ex_seed)  /* Exponential variate generation function. */
{
	/* Return an exponential random variate with mean "mean". */

	return -mean * log(randomvar(ex_seed));
}

float thinning()
{
	float time = time_last_arrival;
	float U1, U2;
	do
		{
			U1 = randomvar(seed1);
			U2 = randomvar(seed2);
			time = time - (1/max_arrival)*log(U1);
			if((int)floor(time / 60) % 2 == 0 && time > hour_coming * 60)
			{
				currArrival = arrival1;
				hour_coming++;
			}
			else if((int)floor(time / 60) % 2 == 1 && time > hour_coming * 60)
			{
				currArrival = arrival2;
				hour_coming++;
			}
			//printf("U1 = %f, U2 = %f, time = %f\n", U1, U2, time);
		}
	while(currArrival/max_arrival < U2);
	time_last_arrival = time;
	return time;
}
float inverse_transform()
{
	float U1 = randomvar(seed1);
	t_star = t_star_last_arrival - log(U1);
	float time;

	//printf ("t_star = %f, lambda = %f, area_lambda=%f, hour_coming=%f, currArrival=%f, sim_time=%f\n", t_star, lambda, area_lambda, hour_coming, currArrival, sim_time);
	if(t_star > threshhold)
	{
		hour_coming++;
		if(currArrival == arrival1)
		{
			threshhold += arrival2*60;
			currArrival = arrival2;
		}
		else
		{
			threshhold += arrival1*60;
			currArrival = arrival1;
		}
	}
	time = (t_star - threshhold + hour_coming*60*currArrival)/currArrival;
	t_star_last_arrival = t_star;
	return time;
}
float randomvar(double *rand_seed)
{      
	hi=floor(*rand_seed/q);
	lo=*rand_seed-q*hi;
	test=a*lo-r*hi;
	if (test>0.0)
		*rand_seed=test;
	else
		*rand_seed=test+m;
	my_random=*rand_seed/m;
	return my_random;
}
